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1.  Introduction 


Coarse-grained  (CG)  potentials  have  the  ability  to  extend  the  length  and  time  scale  that  can  be 
examined  using  simulation.  By  grouping  atoms  together  into  effective  interaction  sites,  it  is 
possible  to  examine  phenomena  and  calculate  macroscopic  properties  inaccessible  from  purely 
fine  resolution  considerations.  Calculations  of  mechanical  properties  in  polymer  systems  are 
particularly  suited  for  this  treatment  because  of  the  range  of  length  and  time  scales  that  describe 
all  the  relevant  physics  within.  For  example,  polymer  bonded  properties  are  of  the  order  of 
0.1  nm,  yet  chains  in  entangled  polymer  systems  may  span  10-100  nm,  both  of  which  affect 
mechanical  properties. 

One  approach  to  investigate  mechanical  properties  of  polymer  systems  has  been  to  use  a  simple 
generic  model  that  incorporates  the  most  pertinent  interactions  within  and  between  atoms.  In 
this  model,  monomers  or  even  groups  of  monomers  are  modeled  as  beads  that  are  connected  to 
adjacent  beads  by  a  finitely  extensible  nonlinear  elastic  (FENE)  potential.  Non-adjacent  beads 
interact  with  a  shifted  Lennard-Jones  potential.  This  potential  was  developed  by  Kremer  and 
Grest  (1)  and  used  to  simulate  strain  hardening  phenomena  (2)  and  glass  transition  (3)  of  linear 
polymers.  The  FENE  potential  has  since  been  modified  slightly  to  allow  for  bond  breaking  and 
has  been  employed  by  many  researchers  to  investigate  a  variety  of  polymer  deformation 
phenomena.  For  example,  Tsige  and  Stevens  studied  adhesion  failure  in  highly  crosslinked 
polymers  (4).  Rottler  et  al.  (5)  and  Panico  et  al.  (6)  applied  tension  on  systems  with  slightly 
different  bonded  potentials  to  study  failure  in  glassy  linear  polymers  at  moderate  and  large 
strains,  respectively.  Mukherji  and  Abrams  also  used  a  modified  FENE  potential  to  examine 
microvoid  formation  and  strain  hardening  in  cross-linked  polymers  (7).  Bulacu  and  van  der 
Giessen  have  changed  polymer  chain  stiffness  by  including  3-  and  4-body  terms  between  beads 
and  studied  its  effect  on  the  glass  transition  (8).  Riggleman  and  coworkers  applied  a  similar 
model  to  examine  mobility  during  creep  of  both  a  linear  polymer  and  nanocomposite,  where 
adjacent  beads  were  connected  with  stiff  harmonic  bonds  instead  of  the  FENE  potential  (9). 
Although  such  studies  have  been  able  to  greatly  advance  our  understanding  of  the  link  between 
microscopic  interactions  and  macroscopic  polymer  deformation  phenomena,  they  lack  the  ability 
to  describe  chemically  specific  systems. 

Many  recipes  exist  to  develop  a  CG  potential  to  model  a  specific  system  of  interest.  In  general, 
chemical  information  can  incorporated  from  atomistic  or  near-atomistic  resolution  simulations  at 
a  given  state  point  to  provide  reference  data.  The  desired  level  of  CG  is  applied  to  this  reference 
data  and  parameterization  is  obtained  by  matching  certain  properties  between  the  two  levels  of 
resolution.  One  procedure  that  deterministically  obtains  CG  potential  parameters  is  called 
iterative  Boltzmann  inversion  (IB I),  which  Reith  and  coworkers  applied  to  parameterize  a  CG 
potential  for  polyisoprene  (JO).  IB  I  sequentially  decreases  the  difference  between  the  free 
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energies  of  the  CG  and  atomistic  distribution  functions.  To  obtain  parameters  for  bonded 
interactions,  one  applies  the  IBI  procedure  using  bond,  angle,  or  dihedral  probability 
distributions,  while  for  nonbonded  interactions  the  radial  distribution  function  is  used.  Miiller- 
Plathe  successfully  parameterized  a  CG  potential  for  polystyrene  at  a  resolution  of  one  CG  bead 
per  monomer  using  IBI  and  examined  the  ability  of  the  potential  to  describe  equilibrium 
properties  at  other  temperatures  (11). 

Another  method  to  obtain  CG  potentials  for  polymer  systems  has  been  used  by  Fritz  and 
coworkers  {12).  Unlike  the  studies  mentioned  previously,  atomistic  simulations  were  not 
performed  in  the  melt.  The  authors  instead  ran  a  constraint  dynamics  simulation  on  two  trimers 
in  vacuum  and  calculated  the  pair  potential  of  mean  force  along  the  distance  coordinate.  The 
potential  of  mean  force  can  then  be  integrated  to  determine  the  pair  potential  between  CG  beads. 
Analogously,  the  bonded  CG  potential  is  calculated  by  obtaining  bond  and  angle  distributions  for 
a  single  long  chain  in  vacuum.  Surprisingly,  the  calculated  potential  was  able  to  reproduce 
densities  and  structural  distributions  in  the  melt  over  a  range  of  temperatures.  Guerrault  and 
coworker  applied  a  similar  method  to  calculate  a  CG  forcefield  for  polyethylene  {13),  yet 
potentials  of  mean  force  were  calculated  in  a  melt  state  from  reference  Isothermal-Isobaric 
(NPT)  Monte  Carlo  simulations.  This  CG  forcefield  was  then  used  to  calculate  equilibrium 
properties  across  a  range  of  chain  lengths. 

Using  the  IBI  approach  or  the  potential  of  mean  force  to  parameterize  a  CG  potential  by 
definition  captures  two  body  correlations,  yet  may  have  severe  limited  state -point  transferability 
{14).  State-point  transferability  refers  to  the  ability  of  a  CG  potential  to  describe  a  given  property 
correctly  at  another  state-point  (i.e.,  temperature,  density,  pressure)  than  the  one  used  in  the 
parameterization  scheme.  Also  of  importance  is  obsen’able  transferability,  which  refers  to  the 
ability  of  a  CG  potential  to  correctly  describe  an  observable  that  has  not  been  used  in  its 
parameterization.  In  addition  to  limited  transferability,  structure -based  approaches  assume  that 
there  is  a  unique  mapping  of  pair  correlations  onto  pairwise  potentials.  The  validity  of  this 
assumption  is  not  readily  apparent  for  complex  systems.  A  direct  method  to  obtain  force  field 
parameters  introduced  by  Izvekov  et  al.  {15)  addresses  the  above  problems.  The  approach 
equates  the  multi-body  force  acting  upon  a  CG  bead  to  that  of  the  atomistically  detailed 
reference  system.  This  bottom-up  multiscale  coarse-graining  (MS-CG)  method  has  been 
successfully  applied  to  ionic  liquids  {16),  small  peptides  {17),  and  molecular  crystals  {18). 

While  the  MS-CG  method  improves  transferability  by  incorporating  multi-body  effects,  no  CG 
method  is  able  to  accurately  reproduce  every  property  of  the  atomistic  ensemble.  Depending 
upon  the  phenomena  of  interest  that  one  needs  to  capture  using  the  CG  potential,  improvement 
can  be  achieved  by  including  additional  information  from  the  atomistic  simulations.  For 
instance,  CG  potentials  parameterized  using  IBI  often  result  in  pressure  or  density  differences 
between  the  CG  and  atomistic  system.  A  common  practice  is  to  shift  the  CG  potential  to  more 
attractive  values  until  pressures  or  densities  match.  Improved  performance  in  equations  of  state 
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calculations  using  MS-CG  has  been  realized  by  aggregating  density  dependent  potentials  to 
investigate  properties  at  vastly  different  state -points  (18,  19). 

While  the  above  methods  to  incorporate  chemical  specificity  have  been  applied  to  calculate  static 
properties,  they  have  been  used  very  infrequently  to  study  mechanical  properties.  Pioneering 
work  conducted  by  Lyulin  and  coworkers  (20)  investigated  the  microscopic  origins  of  strain¬ 
hardening  in  polystyrene  and  polycarbonate.  A  united  atom  model,  where  the  hydrogen  and 
carbon  atoms  are  coarse-grained  into  a  single-site  representation,  was  used.  While  the  united 
atom  model  is  a  very  low  level  of  coarse-graining,  computational  requirements  were  reduced  and 
Lyulin  and  coworkers  were  able  to  demonstrate  that  non-affine  mobility  upon  deformation  is 
linked  to  strain-hardening  effects. 

Finally,  a  recent  paper  by  Majumder  et  al.  (21)  investigated  the  ability  of  extending  the  coarse- 
graining  of  polystyrene  further  and  using  a  coarser  potential  to  calculate  mechanical  properties. 
Like  Lyulin  and  coworkers,  the  authors  used  united  atom  simulations  as  reference  data  upon 
which  the  force-matching  (FM)  technique  of  Izvekov  et  al.  was  applied  to  determine  a  CG 
potential  that  represented  each  polystyrene  monomer  by  one  CG  bead.  The  FM  technique  was 
shown  to  give  a  potential  that  greatly  overestimated  the  pressure.  Consequently,  Majumder  and 
coworkers  modified  the  potential  to  fit  a  Lennard-Jones  type  interaction  and  adjusted  the 
parameters  to  match  the  pressure  between  the  CG  and  united  atom  simulations.  This  procedure 
was  shown  to  reproduce  the  glassy  stress-strain  curve  at  high  deformation  rate. 

In  this  work,  we  further  examine  the  ability  of  chemically  informed  CG  potentials  to  describe 
mechanical  properties.  Molecular  dynamics  simulations  of  amorphous  polymers  for  an  all-atom 
resolution  governed  by  the  PCFF  forcefield  (22)  served  as  reference  data.  The  FM  procedure  of 
Izvekov  et  al.  was  performed  at  a  resolution  of  one  CG  bead  per  monomer.  In  order  to  more 
closely  match  the  structure  between  the  CG  and  atomistic  resolution,  the  iterative  Boltzmann 
inversion  method  was  used  to  determine  the  CG  bonded  potentials.  Stress-strain  extension 
simulations  were  then  performed  and  the  viability  of  the  CG  potentials  assessed. 


2.  Methodology  and  Simulation  Details 


Molecular  dynamics  simulations  were  performed  using  the  open-source  code  Large-scale 
Atomic/Molecular  Massively  Parallel  Simulator  (LAMMPS)  (http://lammps.sandia.gov)  (23). 
The  commercial  visualization  package  Materials  Processes  and  Simulations  (MAPS) 
(http://www.scienomics.com)  was  used  to  create  a  periodic  cell  and  assign  the  PCFF  forcefield 
(22).  Partial  charges  located  on  the  atomic  sites  were  governed  by  the  Condensed -phase 
Optimized  Molecular  Potentials  for  Atomic  Simulation  Studies  (COMPASS)  forcefield,  which 
were  obtained  using  Materials  Studio.  Although  COMPASS  is  an  improvement  upon  the  PCFF 
forcefield,  many  of  the  parameters  are  proprietary  and  cannot  be  ported  to  the  LAMMPS 
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simulation  code.  In  our  molecular  dynamics  simulations  at  the  atomistic  resolution,  we  used  a 

o 

cutoff  of  12.0  A  for  both  Lennard-Jones  (van  der  Waals)  and  Coulombic  nonbonded  interactions. 
Long-range  electrostatic  interactions  were  calculated  using  the  particle-particle  particle-mesh 
method  with  a  precision  of  1.0e-6.  Temperature  was  maintained  by  a  Nose-Hoover  thermostat 
with  a  coupling  constant  of  100  fs.  For  NPT  simulations,  the  pressure  was  controlled  by  a  Nose- 
Hoover  barostat  with  a  coupling  constant  of  1000  fs.  Simulations  at  the  atomistic  resolution  used 
a  timestep  of  1  fs,  while  CG  simulations  used  a  timestep  of  2  fs. 

This  study  focuses  on  one  level  of  CG,  namely  one  CG  site  per  monomer.  To  begin,  an 
atomistic  trajectory  in  both  the  glassy  (298.15  K)  and  rubbery  (500  K)  state  were  generated.  At 
the  outset  of  the  work,  it  was  unclear  which  state-point  would  produce  reference  information  that 
will  best  transfer  to  other  state-points.  One  might  think  that  the  information  obtained  in  the 
glassy  state  would  better  describe  the  physics  necessary  to  obtain  stress-strain  curves  at  other 
glassy  states.  Conversely,  an  atomistic  trajectory  generated  at  500  K  would  explore  more  phase 
space  that  might  produce  a  more  transferrable  forcefield.  To  insure  equilibrated  structures  at  the 
beginning  of  the  glassy  atomistic  trajectory,  initial  amorphous  systems  are  relaxed  and  annealed. 
Relaxation  consists  first  of  energy  minimization  using  the  Polak-Ribiere  conjugate  gradient  with 
a  maximum  of  10000  steps,  and  energy  and  force  tolerance  of  0  and  1.0e-8,  respectively.  After 
minimization,  four  consecutive  microcanonical  ensemble  (NVE)  simulations  of  10000  steps 
were  performed  with  timesteps  of  0.001,  0.01,  0.1,  and  1  fs,  respectively.  Annealing  the 
structure  consisted  of  five  cycles  of  an  NPT  simulation.  Each  cycle  contained  50  ps  at  298.15  K, 
heating  to  600  K  at  a  rate  of  0.67  ps/K,  50  ps  at  600  K,  cooling  to  298.15  K  at  a  rate  of 
0.1656  ps/K,  and  finally  an  energy  minimization.  After  the  annealing  process,  the  structure  was 
subjected  to  10  ns  of  NPT  simulation  under  ambient  conditions,  where  atomic  positions,  forces, 
and  velocities  were  recorded  every  1  ps.  The  rubbery  trajectory  started  from  an  equilibrated 
glassy  system,  which  was  heated  to  500  K  at  a  rate  of  5ps/K  and  then  3  ns  in  the  NPT  ensemble 
were  run,  where  positions,  forces,  and  velocities  were  recorded  every  1  ps. 

A  CG  trajectory  is  mapped  onto  the  last  3  ns  (3000  snapshots)  of  a  trajectory  by  applying  the 
following  equations: 

K,  =  2>,r,  ,  (1) 

i=l,s 

= 1  •  (2) 

i 

where  Ri  is  the  position  of  CG  site  /,  r,  is  the  position  of  atom  i,  and  s  is  the  number  of  atoms  in 
one  CG  site,  c/,  =  where  m,  is  the  mass  of  atom  i  and  MI  is  the  mass  of  CG  site  I.  At  each 
snapshot,  the  atomic  forces  and  velocities  were  coarse-grained  following  an  analogous  equation. 
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(3) 


F,  =  £f, 

i=l,s 

v,  =  5>, 

i=l,s 

Using  a  small  subset  of  the  overall  trajectory  (10  of  the  3000  snapshots),  the  CG  forces  as  a 
function  of  distance  were  calculated  as  a  result  of  matching  to  the  atomistic  data  following  the 
procedure  of  Izvekov  et  al.  (24).  The  resulting  300  functions  for  the  overall  trajectory  were  then 
averaged  together  to  give  the  final  tabulated  CG  force. 

Initially,  it  was  unclear  how  large  the  system  needed  to  be  to  ensure  that  the  properties  calculated 
did  not  suffer  from  system  size  effects.  For  polystyrene,  a  large  system  consisting  of  13  chains 
of  170  monomers  (-35,000  atoms)  and  a  small  system  consisting  of  7  chains  of  20  monomers 
(-2200  atoms)  were  generated.  The  nonbonded  potential  resulting  from  the  force  matching 
procedure  between  beads  for  the  CG  resolution  of  one  site  per  monomer  for  each  system  size  is 
shown  in  figure  1.  The  close  match  between  both  curves  suggests  that  parameterization  of  the 
forcefield  can  be  done  on  a  relatively  small  system,  which  saves  computational  expense. 

Figure  2  displays  the  CG  force  profile  obtained  from  the  glassy  and  rubbery  trajectory  for  a  small 
system  size.  Perhaps  surprisingly,  application  of  the  FM  technique  on  trajectories  generated  at 
two  very  different  temperatures  does  not  produce  large  differences  in  the  potential.  The  largest 
difference  arises  in  smaller  separations,  where  the  curve  parameterized  at  higher  temperature 
corresponds  to  more  repulsive  values.  The  curve  parameterized  at  500  K  also  contains  a  slightly 
deeper  well.  Following  the  procedure  of  Izvekov  et  al.,  we  used  the  500  K  curve  for  the  rest  of 
this  work  because  the  higher  temperature  allows  for  greater  sampling  of  phase  space,  which 
should  allow  for  greater  transferability. 
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Figure  1.  Nonbonded  potential  between  CG  beads  of  polystyrene  parameterized  at  100  K. 
Black  and  blue  curves  correspond  to  results  of  FM  procedure  on 
35,000-  and  2200-atom  systems,  respectively. 


r(  A) 


Figure  2.  Nonbonded  potential  between  CG  beads  of  polystyrene  parameterized  at  two  different 
temperatures.  Blue  and  red  curves  correspond  to  parameterization  at  298  and  500  K, 
respectively. 
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3.  Results  and  Discussion 


To  verify  the  CG  potentials  calculated  through  the  FM  procedure  at  500  K,  bond,  angle,  and 
radial  distribution  functions  of  the  CG  system  were  calculated.  The  CG  system  was  generated 
from  the  last  configuration  of  the  large  atomistic  system  (-35,000  atoms  resulting  in  2210  CG 
sites),  and  distribution  functions  were  calculated  from  a  3-ns  NVT  simulation  with  a  time  step  of 
2  fs  at  the  average  atomistic  density  of  0.88  g/cm3.  Figure  3a-c  displays  the  bond,  angle,  and 
radial  distribution  functions.  The  bonded  and  radial  distributions  match  reasonably  well  for  the 
atomistic  and  CG  resolutions,  while  angular  distributions  are  markedly  different.  Current 
implementation  of  the  FM  technique  is  purely  radial-based,  which  does  not  allow  for  a  direct 
calculation  of  angular  potentials;  thus,  while  1-3  interactions  were  of  the  correct  distance,  their 
orientation  does  not  match  the  atomistic  distributions.  In  addition,  the  potentials  calculated 
exclusively  using  the  FM  technique  result  in  a  density  of  0.78  g/cnr,  an  11%  difference. 


Figure  3.  Distribution  functions  of  the  atomistic  (black)  and  CG  (red)  models,  (a),  (b),  and  (c) 
are  the  bond,  angle,  and  radial  distribution  functions.  Both  bonded  and  nonbonded 
potentials  are  determined  from  a  FM  protocol. 

To  improve  upon  the  internal  structure  of  the  CG  polymer  chains,  the  iterative  Boltzmann 
inversion  method  was  applied  to  calculate  both  bond  and  angle  potentials.  An  initial  guess  of  the 
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potentials  was  based  upon  the  probability  distributions  shown  in  figure  3  a  and  b  (black  curves) 
as  follows: 


U^-kJ  inn,.. 


(4) 


where  i  is  the  bond  distance  or  angle  of  the  distribution.  A  2-ns  simulation  in  the  NVT  ensemble 
(at  the  average  density  from  the  atomistic  trajectory)  was  then  run.  An  improved  potential  was 
then  calculated  from  the  difference  of  the  new  distributions  and  the  target  distributions: 


^ i,new  ^  i, current  ^  ^ 


14  current 
V  ^i, target  J 


(5) 


The  progression  of  the  bond  distribution  upon  the  IBI  procedure  is  shown  in  figure  4. 
Convergence  upon  the  target  distribution  is  efficient,  and  the  CG  bond  distribution  lies  on  top  of 
the  distribution  calculated  from  the  atomistic  distribution  after  three  iterations.  Similar  results 
were  seen  for  the  angular  distribution,  which  along  with  the  bond  and  radial  distribution 
functions,  are  shown  in  figure  5a-c.  The  blend  of  IBI  (for  bonded  interactions)  and  FM  (for 
nonbonded  interactions)  produces  CG  distributions  that  accurately  reproduce  the  atomistic 
structure.  It  should  be  noted  that  the  nonbonded  interactions  between  CG  beads  containing  the 
IBI  bonded  interactions  in  figure  2  is  shifted  to  more  attractive  values  so  that  the  CG  virial  is  5% 
lower  than  the  atomistic  virial.  This  shift  corrects  the  density  discrepancy,  resulting  in  CG 
values  that  are  within  3%  of  the  atomistic  density. 
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Figure  4.  Bond  probability  distribution  for  different  iterations.  Red,  green,  blue  correspond 
to  iteration  1,  2,  and  3,  respectively. 


Figure  5.  Distribution  functions  of  the  atomistic  (black)  and  CG  (red)  models,  (a),  (b),  and 
(c)  are  the  bond,  angle,  and  radial  distribution  functions.  Bonded  and  angular 
potentials  are  determined  using  iterative  Boltzmann  inversion.  The  nonbonded 
potential  is  determined  from  FM  protocol  (same  potential  as  in  figure  2). 
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The  potentials  that  produced  the  distributions  in  figure  5  were  then  used  to  examine  the  CG 
system’s  stress-strain  behavior  in  a  glassy  state.  We  first  examined  the  stress  response  at  100  K 
and  a  strain  rate  of  5e9  s_1  since  these  conditions  were  also  examined  by  Majumder  et  al.  (21). 
While  this  seems  like  an  excessively  low  temperature,  this  ensures  that  the  CG  system  is  indeed 
a  glassy  state.  The  high  strain  rate  should  not  have  any  qualitative  effect  on  the  stress  behavior 
except  to  shift  the  yield  to  higher  values.  Also,  the  purpose  of  this  study  was  not  to  compare  to 
experiment,  but  to  examine  how  well  one  can  describe  stress-strain  behavior  using  CG 
potentials.  To  generate  a  curve  for  the  atomistic  system,  five  independent  systems  were 
constructed.  These  5  independent  systems  of  13  polymers  comprised  of  170  monomers  were 
equilibrated  at  500  K,  quenched  to  100  K  at  a  rate  of  1  ps/K,  and  then  strained  at  a  rate  of  5e9  s_1 
in  the  x-,  y-,  and  z-directions.  Following  the  protocol  of  Vorselaars  and  cowokers  (25),  we 
report  the  stress  as  the  von  Mises  equivalent: 


Every  data  point  is  the  average  of  the  measured  stress  between  each  subsequent  data  point.  In 
addition,  for  the  atomistic  system,  the  data  point  is  the  average  of  each  15  independent  runs.  The 
atomistic  stress-strain  curve  displayed  in  figure  6  has  two  distinct  regimes,  the  elastic  at  small 
strains  and  inelastic  at  large  strains.  The  yield  point  for  the  atomistic  system  at  this  temperature 
and  strain  is  -310  MPa  corresponding  to  -10%  strain  with  a  small  softening  area  corresponding 
to  -30%  strain,  and  finally  a  strain  hardening  regime.  While  this  yield  point  is  very  high,  stress- 
strain  simulations  were  also  run  at  the  conditions  of  Vorselaars  and  coworkers  (25)  study 
(temperature  of  298  K  and  strain  rate  of  le8  s_1),  and  it  was  found  that  our  system’s  mechanical 
properties  matched  very  well  with  theirs.  While  this  validates  the  method  under  which  we 
calculate  the  stress-strain  curves,  the  object  of  this  study  is  to  determine  how  well  a  CG  potential 
can  describe  atomistic  mechanical  properties.  The  middle  curve  in  figure  6  corresponds  to  the 
CG  system  with  bonded  potentials  parameterized  with  IB  I  and  nonbonded  potentials 
parameterized  with  FM,  shifted  slightly  to  reproduce  the  correct  density.  As  evident,  there  is  no 
discernible  yield  point,  along  with  a  substantially  smaller  elastic  and  strain  hardening  modulus. 
The  bottom  curve  corresponds  to  the  CG  system  with  both  bonded  and  nonbonded  potentials 
parameterized  with  FM.  It  should  be  noted  that  the  CG  strain-curves  were  generated  with  only 
one  replica,  yet  adding  the  statistics  from  running  on  the  extra  replicas  would  not  significantly 
change  the  overall  behavior  of  the  stress-strain  curve.  A  larger  CG  system  (160  polymers 
comprised  of  100  monomers)  was  also  simulated  to  ensure  there  are  no  system  size  effects.  The 
results  from  both  CG  systems  were  statistically  equivalent. 
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Figure  6.  Stress-strain  curves  at  100  K  and  a  strain  rate  of  5e9  s  The  top  curve  corresponds 

to  the  fully  atomistic  system  (35,000  atoms).  The  middle  and  bottom  curve  correspond 
to  CG  systems  with  IBI+FM  and  FM-only  potentials,  respectively.  The  CG  systems 
are  the  same  absolute  sizes  as  the  atomistic  system  (2200  CG  sites). 


It  is  apparent  from  figure  6  that  the  pertinent  physics  is  missing  from  the  CG  system  to 
accurately  describe  stress -strain  behavior  found  in  the  atomistic  system.  It  is  well  known  that 
CG  potentials  result  in  faster  diffusion  times  of  CG  particles  as  compared  to  atomistic 
counterparts  because  they  are  smoother  and  softer  (26).  One  way  of  controlling  the  diffusion  of 
particles  in  a  simulation  is  to  run  with  a  dissipative  particle  dynamics  (DPD)  thermostat.  In 
addition  to  the  force  arising  from  dispersion  and/or  electrostatic  interactions,  the  DPD  thermostat 
includes  dissipative  and  random  forces.  The  arithmetic  forms  of  these  two  forces  are  given 
below: 


K 
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-r  1 
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°°ij 


Jdt 


r  r.y 

\-^~ 
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C  J 


(7) 


where  y,  <j,  6tj,  and  St  are  the  dissipative  scaling  factor,  noise  level,  random  number,  and 
timestep,  respectively.  Espanol  and  Warren  derived  the  required  fluctuation-dissipation  relation 
between  the  dissipative  scaling  factor  and  noise  level  to  sample  from  the  canonical  ensemble 
(27).  By  setting  the  dissipative  term,  one  can  effectively  slow  down  the  system  dynamics,  and  in 
terms  of  CG  systems,  match  the  dynamics  of  the  atomistic  system. 
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In  a  glassy  system  it  is  very  difficult  to  measure  diffusion  because  it  is  exceedingly  small.  We 
therefore  measured  the  diffusion  at  higher  temperature  and  then  extrapolated  to  lower 
temperatures.  In  addition,  we  believe  that  stress-strain  phenomena  captured  in  our  atomistic 
systems  are  dependent  mostly  on  bead-bead  interactions,  and  not  on  entanglement.  Therefore  we 
are  interested  in  matching  the  diffusion  on  the  bead  length  scale  and  not  on  the  polymer  length 
scale.  To  match  diffusion  on  this  scale,  we  simulate  a  system  of  500  molecules  consisting  of 
3  monomers  of  polystyrene.  The  trimer  contains  all  of  the  interaction  components  of  a  longer 
polymer  chain,  i.e.,  bond,  angle,  and  nonbond.  Figure  7  displays  the  dependence  of  the  CG 
trimer  diffusion  with  respect  to  dissipative  scaling  factor  y  at  400  K.  The  nonbonded  potential 
used  in  the  CG  simulation  was  parameterized  using  the  FM  method  with  the  atomistic  trimer 
diffusive  trajectory  at  400  K  as  input,  while  the  IBI  method  was  used  to  obtain  the  angle  and 
bond  potentials.  We  have  fit  an  inverse  power  law  to  the  data  points  to  estimate  the  y  value  at 
which  the  CG  and  atomistic  diffusion  coefficient  match.  This  procedure  was  repeated  at  higher 
temperatures  with  the  results  presented  in  figure  8. 


Figure  7.  Relative  diffusive  speedup  of  the  CG  trimer  system  at  400  K  as  a  function  of  the 

dissipative  term  y.  y  is  in  units  of  kcal  /(mol  fs).  The  symbols  are  simulation  results, 
while  the  red  curve  is  a  fit  to  an  inverse  power  law.  An  arrow  indicates  the  value  at 
which  the  CG  and  atomistic  systems  have  the  same  diffusion. 
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Figure  8.  Dissipative  term  rat  which  the  diffusion  of  atomistic  and  CG  systems  are  equal. 

The  symbols  are  simulation  results,  while  the  red  curve  is  a  fit  to  a  cubic  polynomial. 

A  cubic  polynomial  was  fit  to  the  data  and  this  equation  was  used  to  extrapolate  to  a  value  of 
y=  45500  kcal/(mol  fs)  at  100  K.  Performing  stress-strain  simulations  with  this  dissipative  term 
results  in  a  curve  presented  in  figure  9,  which  greatly  improves  the  CG  models’s  ability  to  match 
the  stress-strain  behavior  compared  to  the  atomistic  system.  The  elastic  modulus  between  the 
CG  model  determined  using  DPD  and  the  atomistic  model  are  now  equal  within  statistical 
accuracy;  however,  the  inelastic  region  does  not  match  as  well.  Moreover,  the  CG  system  does 
not  show  any  yield  point  and  exhibits  a  strain  hardening  modulus  that  is  larger  than  the  atomistic 
system. 

While  figure  9  displays  a  systematic  parameterization  of  a  CG  forcefield  to  describe  mechanical 
properties,  the  potentials  must  have  some  transferability  to  other  strain  rates  and  temperatures  to 
be  practical.  To  probe  the  robustness  of  the  parameters  used  in  the  DPD  thermostat,  we  ran 
simulations  at  a  slightly  slower  strain  rate  of  le9  s_1  and  100  K,  as  well  as  the  same  strain  rate 
and  a  higher  temperature  of  298  K,  which  are  presented  in  figure  10.  For  the  slightly  slower 
strain  rate,  the  dissipative  term  remained  constant,  while  at  298  K  y=  8800  kcal/(mol  fs).  Strains 
were  performed  on  only  one  replica;  thus,  the  data  in  figure  10  is  a  little  noisier  than  the  other 
stress-strain  curves. 
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Figure  9.  Stress-strain  curves  at  100  K  and  a  strain  rate  of  5e9  s  .  The  black  and  red  curves 

are  the  same  as  in  figure  6,  while  the  blue  curve  was  generated  using  a  DPD  thermostat 
with  a  dissipative  term  y=  45500  kcal/(mol  fs). 

For  the  slightly  slower  strain  rate  of  le9  s  ,  the  atomistic  stress-strain  behavior  is  essentially 
identical  to  the  higher  strain  rate,  yet  the  CG  system  displays  marked  differences  between  the 
two  strain  rates.  At  the  slower  strain  rate,  the  CG  system  has  a  much  lower  elastic  modulus  and 
the  behavior  is  almost  rubbery  in  nature.  The  atomistic  stress-strain  curve  at  298  K  still  appears 
glassy  in  nature  with  a  smaller  elastic  modulus  than  at  100  K,  as  is  expected.  The  CG  stress- 
strain  curve  at  the  higher  temperature  is  very  similar  to  the  higher  strain  rate-lower  temperature 
curve,  with  no  glassy  characteristics.  Overall,  figure  10  suggests  that  the  friction  coefficient 
necessary  to  reconcile  the  CG  and  atomistic  mechanical  properties  may  be  both  temperature  and 
strain-rate  dependent. 
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Figure  10.  Stress-strain  curves  at  two  different  strain  rates  and  temperatures.  Solid  lines 

correspond  to  atomistic  systems  and  dashed  lines  correspond  to  CG  systems.  Black  curves 
are  at  a  strain  rate  of  5e9  s  and  100  K  (for  CG  system  y=  45500  kcal/(mol  fs)).  Blue 
curves  are  at  a  strain  rate  of  le9  s  1  and  100  K  (for  CG  system  y=  45500  kcal/(mol  fs)). 

Red  curves  are  at  a  strain  rate  of  5e9  s  and  298  K  (for  CG  system  y=  8800  kcal/(mol  fs)). 

While  adding  dissipative  forces  to  the  system  via  the  thermostat  is  not  amenable  to  simple 
parameterization,  an  alternative  is  to  directly  alter  the  nonbonded  interactions.  Because  the 
modulus  of  the  CG  system  is  much  smaller  than  the  atomistic  system,  it  is  apparent  that  the 
nonbonded  potential  needs  to  be  more  attractive  to  obtain  a  stiffer  response.  The  simplest 
approach  to  induce  more  attraction  is  to  uniformly  shift  the  nonbonded  forces  to  lower  values. 
Following  this  approach,  various  constants  were  added  to  the  nonbonded  interaction  and  for  each 
value  the  system  was  equilibrated  at  500  K,  quenched  at  1  K/ps,  and  strained  to  10%  in  the  x- 
direction.  For  the  value  that  demonstrated  an  elastic  modulus  closest  to  the  atomistic  model,  the 
full  stress-strain  curve  was  generated  in  each  direction.  Figure  1 1  displays  both  the  stress-strain 
curves  corresponding  to  the  system  with  the  shifted  nonbonded  potential  that  matches  the 
atomistic  elastic  modulus  (red  curve),  as  well  a  system  with  a  slightly  larger  shift  (blue  curve). 
The  shifted  potential  that  matches  the  atomistic  system  contains  a  minimum  of  -0.93  kcal/mol  at 

o  o 

7.866  A,  down  from  -0.53  kcal/mol  at  8.1775  A  (see  inset  to  figure  11).  The  overly  stiff  system 

o 

corresponds  to  a  potential  containing  a  minimum  of  -1.04  kcal/mol  at  7.86  A. 
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Figure  11.  Stress-strain  curves  at  100  K  and  a  strain  rate  of  5e9  s  The  black  curve  corresponds 
to  the  atomistic  system.  The  green,  red,  and  blue  curves  are  for  CG  systems  containing 
nonbonded  interactions  displayed  in  the  inset.  The  green  curve  is  the  same  as  the  lower 
curve  in  figure  9. 

While  the  elastic  region  below  5%  strain  can  be  matched  by  a  constant  shift  of  the  force,  the  CG 
potential  is  unable  to  reproduce  the  atomistic  yield  point,  with  a  deviation  occurring  at  a  stress  of 
-200  MPa.  In  addition,  the  shifted  potential  significantly  alters  the  structure  of  the  CG  system. 
Figure  12  displays  the  resulting  bond,  angle,  and  radial  distributions  upon  equilibration  at  500  K. 
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Figure  12.  Distribution  functions  of  the  atomistic  (black)  and  CG  (red)  models,  (a),  (b),  and  (c) 
are  the  bond,  angle,  and  radial  distribution  functions. 


The  distributions  in  figure  12  demonstrate  the  destruction  of  the  structure  resulting  from  the 
increase  in  attraction  between  CG  beads.  The  shifted  nonbonded  potential  results  in  chains 
adopting  smaller  interbead  separations,  as  well  as  more  acute  angular  orientations.  The  radial 
distribution  function  also  displays  a  much  sharper  peak  at  smaller  separations  as  compared  to  the 
atomistic  system. 


While  shifting  the  nonbonded  potential  to  more  attractive  values  has  the  desired  effect  on  the 
mechanical  properties,  it  has  a  decidedly  undesired  effect  on  the  structure  of  the  CG  system. 
From  the  inset  in  figure  1 1,  it  is  apparent  that  shifting  the  force  by  a  constant  amount  influences 
the  potential  over  the  entire  range  of  interactions.  Alternatively,  we  attempted  to  lessen  the 
destruction  of  the  CG  structure  by  shifting  the  force  in  only  a  small  subset  of  the  overall 
interaction  range.  We  thus  kept  the  location  of  the  zero  and  minimum  in  the  nonbonded 
potential  constant  by  adding  a  sinusoidal  function  to  the  force 
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where  ro  is  the  radius  at  which  the  nonbonded  potential  equals  0,  A  is  the  magnitude  of  the  shift 
to  the  nonbonded  potential,  and  w  is  twice  the  distance  between  the  0  and  minimum  in  the 
nonbonded  potential.  The  only  adjustable  parameter  in  the  above  equations  is  A,  which  controls 
the  depth  of  the  potential  energy  minimum,  but  not  the  location.  Analogous  to  the  procedure 
above,  various  values  of  A  were  chosen  and  the  new  potentials  calculated.  For  each  new 
potential,  a  CG  system  was  equilibrated  at  500  K  for  2  ns,  quenched  to  100  K  at  1  K/ps,  and  then 
strained  to  10%.  Figure  13  displays  the  full  stress-strain  curve  corresponding  to  the  choice  of  A 
that  most  satisfactorily  matches  the  atomistic  curve  the  best.  The  inset  contains  the  nonbonded 
potential  that  corresponds  to  shifting  the  force  with  the  sinusoidal  function  of  equation  8. 
Compared  to  the  potential  obtained  by  shifting  the  force  by  a  constant  value,  the  sinusoidal  shift 
localizes  the  increase  in  attraction  to  a  much  smaller  range.  The  locally  shifted  potential  results 
in  a  yield  point  and  strain  hardening  modulus  that  more  closely  matches  the  atomistic  system  as 
compared  to  the  uniformly  shifted  potential  found  in  figure  11.  Figure  14  displays  the  structural 
distribution  functions  associated  with  this  locally  shifted  potential. 


Figure  13.  Stress-strain  curves  at  100  K  and  a  strain  rate  of  5e9  s  The  black  curve  corresponds 

to  the  atomistic  system.  The  green,  red,  and  purple  curves  are  for  CG  systems  containing 
nonbonded  interactions  displayed  in  inset.  The  green  and  red  curves  are  the  same  as  the 
lower  curve  in  figure  1 1 . 
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Figure  14.  Distribution  functions  of  the  atomistic  (black)  and  CG  (red)  models,  (a),  (b),  and  (c) 
are  the  bond,  angle,  and  radial  distribution  functions. 

While  the  locally  shifted  potentials  produce  structural  distributions  that  deviate  from  the 
atomistic  distributions,  when  compared  to  the  uniformly  shifted  potentials,  the  deviation  is  not  as 
significant.  The  CG  bonded  and  angle  distributions  are  shifted  to  lower  values  because  of  the 
added  attraction,  yet  the  radial  distribution  is  affected  only  locally.  A  peak  appears  at  slightly 
larger  separations  than  the  atomistic  distribution.  Overall,  the  protocol  to  add  attraction  to  the 
CG  model  locally  rather  than  uniformly  results  in  slightly  better  stress -strain  behavior  and 
significantly  better  structural  results. 

The  locally  shifted  potential  was  able  to  minimally  impact  equilibrium  structural  properties  while 
matching  mechanical  properties  at  100  K  and  a  high  strain  rate,  however,  we  also  examined  its 
applicability  to  other  strain  rates  and  temperatures.  Figure  15  displays  stress-strain  curves  at  a 
slower  strain  rate  of  le8  s_1  and  a  temperature  of  100  K.  The  locally  shifted  potential  matches 
the  atomistic  behavior  in  the  elastic  region,  yet  yields  at  a  slightly  smaller  value.  In  the  elastic 
regime  the  locally  shifted  potential  and  atomistic  systems  display  very  similar  behavior  as  well. 
The  uniformly  shifted  potential  yields  at  a  much  smaller  value  than  the  locally  shifted  CG  system 
and  atomistic  system,  and  appears  to  have  two  distinct  strain-hardening  regions.  From  0.05  to 
~0.4  strain,  the  uniformly  shifted  potential  has  a  region  with  a  large  strain  hardening  modulus, 
while  strains  larger  than  0.4  produces  a  modulus  very  similar  in  value  to  the  atomistic  system. 
The  fluctuations  in  the  stress  of  the  uniformly  shifted  potential  are  much  larger  in  the  elastic 
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region  compared  to  the  locally  shifted  potential  system.  Note  that  the  CG  stress-strain  curves 
were  generated  using  only  one  replica,  thus  it  is  not  appropriate  to  compare  the  CG  system 
fluctuations  to  the  atomistic  system. 


Figure  15.  Stress-strain  curves  at  100  K  and  a  strain  rate  of  le8  s  The  black  curve  corresponds 

to  the  atomistic  system.  The  red  (uniformly  shifted  potential)  and  purple  (locally  shifted 
potential)  curves  are  for  CG  systems  containing  nonbonded  interactions  displayed  in  the 
inset  of  figure  12. 

Finally,  we  examined  the  transferability  of  the  shifted  potentials  to  describe  stress-strain 
behavior  at  a  different  temperature.  Figure  16  displays  the  stress-strain  curves  at  298  K  and 
5e9  s  .  The  locally  shifted  potential  system  has  a  smaller  modulus  and  yield  point  compared  to 
the  atomistic  system.  While  the  elastic  region  of  the  locally  shifted  potential  CG  system  does  not 
match  the  atomistic  system,  the  strain-hardening  behavior  does  satisfactorily  agree.  The 
uniformly  shifted  potential  transfers  better  to  the  higher  temperature  than  the  locally  shifted 
potential,  with  both  the  elastic  and  inelastic  region  reasonably  matching  the  atomistic  model 
behavior. 
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Figure  16.  Stress-strain  curves  at  298  K  and  a  strain  rate  of  5e9  s  The  black  curve  corresponds 

to  the  atomistic  system.  The  red  (uniformly  shifted  potential)  and  purple  (locally  shifted 
potential)  curves  are  for  CG  systems  containing  nonbonded  interactions  displayed  in  the 
inset  of  figure  12. 


4.  Conclusion 


We  investigated  the  ability  of  a  chemically  informed  CG  potential  of  a  single  resolution  (one 
monomer  to  one  CG  particle)  to  describe  the  mechanical  behavior  at  glassy  temperatures  and 
high  strain  rate.  The  FM  technique  of  Izvekov  et  al.  was  used  to  determine  nonbonded 
interactions,  while  the  IBI  method  was  used  to  determine  bonded  and  angular  interactions. 
Applied  in  a  straightforward  manner,  these  techniques  were  insufficient  to  reproduce  stress-strain 
behavior.  Slowing  down  the  dynamics  of  the  CG  system  by  implementing  the  constant 
temperature  DPD  method  proved  intractable.  Matching  the  diffusion  coefficients  at  higher 
temperatures  by  altering  the  DPD  friction  parameter  y  gave  satisfactory  mechanical  data  at  100  K 
and  a  strain  rate  of  5e9  s_1,  yet  did  not  transfer  to  different  strain  rates  and  temperatures.  Two 
methods  to  introduce  more  attractive  nonbonded  interactions  to  the  CG  system  were  explored,  a 
uniformly  shifted  force  and  a  locally  shifted  force.  The  locally  shifted  force  (which  retained  the 
potential’s  zero  and  minimum  radial  location)  produced  a  better  yield  point  as  well  as  smaller 
fluctuations  in  the  inelastic  region  as  compared  to  the  uniformly  shifted  force.  In  addition,  the 
locally  shifted  potential  has  minimal  impact  on  the  equilibrium  structure  of  the  system.  While 
transferability  to  different  strain  rates  was  reasonable,  stress-strain  behavior  at  higher 
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temperatures  was  poor  compared  to  the  uniformly  shifted  force.  Future  studies  will  examine  a 
higher  resolution  CG  model;  two  CG  particles  per  one  polymer  monomer.  The  objective  will  be 
to  determine  the  quantity  (if  any)  that  the  potential  needs  to  be  altered  to  recover  the  atomistic 
models  mechanical  behavior.  In  addition,  dihedral  potentials  will  be  added  to  this  higher 
resolution  and  its  effect  will  be  explored. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


CG 

COMPASS 

DPD 
FENE 
FM 
IB  I 

LAMMPS 

MAPS 

MS-CG 

NPT 

NVE 


coarse-grained 

Condensed-phase  Optimized  Molecular  Potentials  for  Atomic  Simulation 
Studies 

dissipative  particle  dynamics 
finitely  extensible  nonlinear  elastic 
force-matching 
iterative  Boltzmann  inversion 

Large-scale  Atomic/Molecular  Massively  Parallel  Simulator 

Materials  Processes  and  Simulations 

multiscale  coarse-graining 

Isothermal-Isobaric 

microcanonical  ensemble 
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